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Open Source License 


Copyright @ 2015, Illinois Institute of Technology. All rights reserved. 

Redistribution and use in source and binary forms, with or without modification, are permitted provided 
that the following conditions are met: 


• Redistributions of source code must retain the above copyright notice, this list of conditions and the 
following disclaimer. 


• Redistributions in binary form must reproduce the above copyright notice, this list of conditions and 
the following disclaimer in the documentation and/or other materials provided with the distribution. 

• Neither the name of Illinois Institute of Technology nor the names of its contributors may be used to 
endorse or promote products derived from this software without specihc prior written permission. 

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER AND CONTRIBUTORS “AS IS” 
AND WITHOUT ANY WARRANTY OF ANY KIND, WHETHER EXPRESS, IMPLIED, STATUTORY 
OR OTHERWISE, INCLUDING WITHOUT LIMITATION WARRANTIES OF MERCHANTABILITY, 
FITNESS FOR A PARTICULAR USE AND NON-INERINGEMENT, ALL OF WHICH ARE HEREBY 
EXPRESSLY DISCLAIMED. MOREOVER, THE USER OF THE SOETWARE UNDERSTANDS AND 
AGREES THAT THE SOFTWARE MAY CONTAIN BUGS, DEFECTS, ERRORS AND OTHER PROB¬ 
LEMS THAT COULD CAUSE SYSTEM FAILURES, AND ANY USE OF THE SOETWARE SHALL BE 
AT USER?S OWN RISK. THE COPYRIGHT HOLDERS AND CONTRIBUTORS MAKES NO REPRE¬ 
SENTATION THAT THEY WILL ISSUE UPDATES OR ENHANCEMENTS TO THE SOFTWARE. 

IN NO EVENT WILL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE EOR ANY DI¬ 
RECT, INDIRECT, SPECIAL, INCIDENTAL, CONSEQUENTIAL, EXEMPLARY OR PUNITIVE DAM¬ 
AGES, INCLUDING, BUT NOT LIMITED TO, DAMAGES FOR INTERRUPTION OF USE OR FOR 
LOSS OR INAGGURAGY OR CORRUPTION OF DATA, LOST PROFITS, OR COSTS OF PROCURE¬ 
MENT OE SUBSTITUTE GOODS OR SERVICES, HOWEVER CAUSED (INCLUDING BUT NOT LIM¬ 
ITED TO USE, MISUSE, INABILITY TO USE, OR INTERRUPTED USE) AND UNDER ANY THEORY 
OE LIABILITY, INGLUDING BUT NOT LIMITED TO GONTRAGT, STRIGT LIABILITY, OR TORT 
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OE THE USE OE THIS 
SOETWARE, EVEN IE ADVISED OF THE POSSIBILITY OF SUGH DAMAGE AND WHETHER OR 
NOT THE GOPYRIGHT HOLDER AND CONTRIBUTORS WAS OR SHOULD HAVE BEEN AWARE 
OR ADVISED OF THE POSSIBILITY OE SUCH DAMAGE OR EOR ANY GLAIM ALLEGING INJURY 
RESULTING EROM ERRORS, OMISSIONS, OR OTHER INAGGURAGIES IN THE SOFTWARE OR 
DESTRUGTIVE PROPERTIES OF THE SOFTWARE. TO THE EXTENT THAT THE LAWS OE ANY 
JURISDIGTIONS DO NOT ALLOW THE FOREGOING EXGLUSIONS AND LIMITATION, THE USER 
OE THE SOETWARE AGREES THAT DAMAGES MAY BE DIEEIGULT, IF NOT IMPOSSIBLE TO 
CALCULATE, AND AS A RESULT, SAID USER HAS AGREED THAT THE MAXIMUM LIABILITY 
OE THE GOPYRITGHT HOLDER AND GONTRIBUTORS SHALL NOT EXGEED US$100.00. 

THE USER OE THE SOFTWARE ACKNOWLEDGES THAT THE SOETWARE IS BEING PROVIDED 
WITHOUT GHARGE, AND AS A RESULT, THE USER, ACKNOWLEDGING THAT HE OR SHE HAS 
READ THE SAME, AGREES THAT THE FOREGOING LIMITATIONS AND RESTRIGTIONS REP¬ 
RESENT A REASONABLE ALLOGATION OE RISK. 
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1 Introduction 


Automatic and adaptive approximation, optimization, or integration of functions in a cone with guarantee 
of accuracy is a relatively new paradigm [7]. Our purpose is to create an open-source MATLAB package, 
Guaranteed Automatic Integration Library (GAIL) [5], following the philosophy of reproducible research 
championed by Claerbout [6] and Donoho [1] , and sustainable practices of robust scientific software develop¬ 
ment [12]. For our conviction that true scholarship in computational sciences are characterized by reliable 
reproducibility [3, 4, 2], we employ the best practices in mathematical research and software engineering 
known to us and available in MATLAB. 

The rest of this document describes the key features of functions in GAIL, which includes one-dimensional 
function approximation [7, 8] and minimization [14] using linear splines, one-dimensional numerical integra¬ 
tion using trapezoidal rule [7], and last but not least, mean estimation and multidimensional integration by 
Monte Carlo methods [9, 11] or Quasi Monte Carlo methods [13, 10]. 


1.1 Downloads 

GAIL can be downloaded from 
http;//code.google.com/p/gail/ 

Alternatively, you can get a local copy of the GAIL repository with this command: 
git clone https://github.com/GailGithub/GAlL_Dev.git 

1.2 Requirements 

You will need to install MATLAB 7 or a later version. 

1.3 Documentation 

Detailed documentation is available at GAIL_Matlab/Documentation. 

1.4 General Usage Notes 

GAIL Version 2.1 [5] includes the following eight algorithms: 

1. funappx_g [7, 8]: One-dimensional function approximation on bounded interval 

2. funmin_g [14]: global minimum value of univariate function on a closed interval 

3. integraLg [7]: One-dimensional integration on bounded interval 

4. meanMC_g [9]: Monte Carlo method for estimating mean of a random variable 

5. meanMCBer_g [11]: Monte Carlo method to estimate the mean of a Bernoulli random variable 

6. cubMC_g [9]: Monte Carlo method for numerical multiple integration 

7. cubLattice_g [13]: Quasi-Monte Carlo method using rank-1 Lattices cubature for a d-dimensional 
integration 

8. cubSoboLg [10]: Quasi-Monte Carlo method using Sobol’ cubature for a d-dimensional integration 
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1.5 Installation Instruction 


1. Unzip the contents of the zip file to a directory and maintain the existing directory and subdirectory 
structure. (Please note: If you install into the toolbox subdirectory of the MATLAB program hierarchy, you 
will need to click the button “Update toolbox path cache” from the File/Preferences... dialog in MATLAB.) 


2. In MATLAB, add the GAIL directory to your path. This can be done by running GAIL_Install .m. 
Alternatively, this can be done by selecting “File/Set Path...” from the main or Command window menus, 
or with the command pathtool. We recommend that you select the “Save” button on this dialog so that 
GAIL is on the path automatically in future MATLAB sessions. 

3. To check if you have installed GAIL successfully, type help funappx_g to see if its documentation shows 
up. 

Alternatively, you could do this: 

1. Download DownloadInstallGaiL2_l.m and put it where you want GAIL to be installed. 

2. Execute it in MATLAB. 

To uninstall GAIL, execute GAlL_Uninstall. 

To reinstall GAIL, execute GAlL_lnstall. 


1.6 Tests 

We provide quick doctests for each of the functions above. To run doctests in funappx_g, for example, issue 
the command doctest funappx^. 

We also provide unit tests for MATLAB version 8 or later. To run unit tests for funmin^, for instance, 
execute run(ut Junmin_g). 


1.7 Contact Information 

Please send any queries, questions, or comments to 
gail-users@googlegroups.com 


1.8 Website 

For more information about GAIL, visit GAIL Project website. 
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2 funappx_g 

l-D guaranteed locally adaptive function approximation (or function recovery) on [a,b] 

2.1 Syntax 

fappx = funappx_g(f) 

fappx = funappx_g(f,a,b,abstol) 

fappx = funappx_g(f,’a’,a,’b’,b,’abstor,abstol) 

fappx = funappx_g(f,in_param) 

[fappx, out.param] = funappx_g(f,...) 

2.2 Description 

fappx = funappx_g(f) approximates function f on the default interval [0,1] by an approximated function 
handle fappx within the guaranteed absolute error tolerance of le-6. When Matlab version is higher or equal 
to 8.3, fappx is an interpolant generated by griddedinterpolant. When Matlab version is lower than 8.3, 
fappx is a function handle generated by ppval and interpl. Input f is a function handle. The statement y = 
f(x) should accept a vector argument x and return a vector y of function values that is of the same size as 

X. 

fappx = funappx_g(f,a,b,abstol) for a given function f and the ordered input parameters that define the 
finite interval [a,b], and a guaranteed absolute error tolerance abstol. 

fappx = funappx_g(f,’a’,a,’b’,b,’abstol’,abstol) approximates function f on the finite interval [a,b], given a 
guaranteed absolute error tolerance abstol. All four field-value pairs are optional and can be supplied in 
different order. 

fappx = funappx_g(f,in_param) approximates function f on the finite interval [in_param.a,in_param.b], given 
a guaranteed absolute error tolerance in_param.abstol. If a field is not specified, the default value is used. 

[fappx, out_param] = funappx_g(f,...) returns an approximated function fappx and an output structure 
out_param. 

Input Arguments 

• f — input function 

• in_param.a — left end point of interval, default value is 0 

• in_param.b — right end point of interval, default value is I 

• in_param.abstol — guaranteed absolute error tolerance, default value is le-6 

Optional Input Arguments 

• in_param.nlo — lower bound of initial number of points we used, default value is 10 

• in_param.nhi — upper bound of initial number of points we used, default value is 1000 

• in_param.nmax — when number of points hits the value, iteration will stop, default value is le7 
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• in_param.maxiter — max number of iterations, default value is 1000 


Output Arguments 


• fappx — approximated function handle (Note: When Matlab version is higher or equal to 8.3, fappx 
is an interpolant generated by griddedinterpolant. When Matlab version is lower than 8.3, fappx is a 
function handle generated by ppval and interpl.) 

• out_param.f — input function 

• out_param.a — left end point of interval 

• out_param.b — right end point of interval 

• out_param.abstol — guaranteed absolute error tolerance 

• out_param.nlo — a lower bound of initial number of points we use 

• out_param.nhi — an upper bound of initial number of points we use 

• out_param.nmax — when number of points hits the value, iteration will stop 

• out_param.maxiter — max number of iterations 

• out_param.ninit — initial number of points we use for each sub interval 

• out_param.exit — this is a number defining the conditions of success or failure satisfied when finishing 
the algorithm. The algorithm is considered successful (with out_param.exit == 0) if no other flags 
arise warning that the results are certainly not guaranteed. The initial value is 0 and the final value 
of this parameter is encoded as follows: 1 If reaching overbudget. It states whether the max budget 
is attained without reaching the guaranteed error tolerance. 

2 If reaching overiteration. It states whether the max iterations is attained without reaching the 
guaranteed error tolerance. 

• out_param.iter — number of iterations 

• out_param.npoints — number of points we need to reach the guaranteed absolute error tolerance 

• out_param.errest — an estimation of the absolute error for the approximation 

• out_param.nstar — final value of the parameter defining the cone of functions for which this algorithm 
is guaranteed for each subinterval; nstar = ninit-2 initially 

2.3 Guarantee 

For [a, 6], there exists a partition 

P = {[io,ii], [^ 1 ,^ 2 ],- • ■, [tL-i,tL]},a = to <ti < ■ ■ ■ < tL = b. 

If the function to be approximated, / satisfies the cone condition 

ti-ti^i 00 

for each sub interval where 1 < I < L, then the fappx [output by this algorithm is guaranteed to 

satisfy 

11/ - fappx\\oo < abstol. 


< 


2nstar 
ti — ti-i 
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2.4 Examples 

Example 1 

f = @(x) x.~2; [fappx, out_parain] = funappx_g(f) 

7o Approximate function x~2 with default input parameter to make the error 
7„ less than le-6. For MATLAB version 8.3 onwards, we see: 


fappx = 


griddedlnterpolant with properties: 


GridVectors: 
Values: 
Method: 
ExtrapolationMethod: 


{[1x3169 double]} 
[1x3169 double] 
'linear' 

'linear' 


out_param = 


f: 

(§(x)x.‘2 

a: 

0 

b: 

1 

abstol: 

l.OOOOe-06 

nlo: 

10 

nhi: 

1000 

nmax: 

10000000 

maxiter: 

1000 

ninit: 

100 

exit: 

[2x1 logical] 

iter: 

6 

npoints: 

3169 

errest: 

2.7429e-07 

nstar: 

[1x32 double] 


For earlier versions of MATLAB, we have: 


fappx = 

@(x)ppval(pp,x) 


.param = 

f: 

(§(x)x.‘2 

a: 

0 

b: 

1 

abstol: 

l.OOOOe-06 

nlo: 

10 

nhi: 

1000 

nmax: 

10000000 
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maxiter: 
ninit: 
exit: 
iter: 
npoints: 
errest; 
nstar: 


1000 

100 

[2x1 logical] 

6 

3169 

2.7429e-07 

[10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 . . . 
10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 ] 


Example 2 

[fappx, out_parajn] = funappx_g(@(x) x.~2,0,100,le-7,10,1000,leS) 

7» Approximate function x"2 on [0,100] with error tolerance le-7, cost 
7, budget 10000000, lower bound of initial number of points 10 and upper 
7, bound of initial number of points 100 

fappx = 


griddedlnterpolant with properties: 


GridVectors: 
Values: 
Method: 
ExtrapolationMethod: 


■[[1x977921 double]} 
[1x977921 double] 

’linear’ 

’linear’ 


out_param = 


a: 

0 

abstol: 

l.OOOOe-07 

b: 

100 

f: 

(§(x)x.‘2 

maxiter: 

1000 

nhi: 

1000 

nlo: 

10 

nmax: 

100000000 

ninit: 

956 

exit: 

[2x1 logical] 

iter: 

11 

npoints: 

977921 

errest: 

3.7104e-08 

nstar: 

[1x1024 double] 


Example 3 

clear in_parcun; in_param.a = -20; in_parcmi.b = 20; in_param.nlo = 10; 
in_param.nhi = 100; in_param.nmax = leS; in_param.abstol = le-7; 
[fappx, out_param] = funappx_g(@(x) x.~2, in_param) 

7o Approximate function x"2 on [-20,20] with error tolerance le-7, cost 
7o budget 1000000, lower bound of initial number of points 10 and upper 
7o bound of initial number of points 100 
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fappx = 


griddedinterpolant with properties: 


GridVectors: 
Values: 
Method: 
ExtrapolationMethod: 


{[1x385025 double]} 
[1x385025 double] 
’linear’ 

’linear’ 


out_parain = 


a: 

-20 

abstol: 

l.OOOOe-07 

b: 

20 

f: 

@(x)x.‘2 

maxiter: 

1000 

nhi: 

100 

nlo: 

10 

nmax: 

100000000 

ninit: 

95 

exit: 

[2x1 logical] 

iter: 

13 

npoints: 

385025 

errest: 

2.6570e-08 

nstar: 

[1x4096 double] 


Example 4 


clear in_parain; f = @(x) x."2; 

[fappx, out_parain] = funappx_g(f,’a’,-10,’b’,50,’nmax’,le6,’abstol’,le-7) 

7o Approximate function x~2 with error tolerance le-7, cost budget 1000000, 
7„ lower bound of initial number of points 10 and upper 
7„ bound of initial number of points 100 


fappx = 


griddedinterpolant with properties: 


GridVectors: 
Values: 
Method: 
ExtrapolationMethod: 


{[1x474625 double]} 
[1x474625 double] 

’linear’ 

’linear’ 


out_param = 


a: -10 

abstol: l.OOOOe-07 
b: 50 

f: (a(x)x.‘2 


maxiter: 1000 
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nhi: 
nlo: 
nmax: 
ninit: 
exit: 
iter: 
npoints; 
errest: 
nstar: 


1000 

10 

1000000 

928 

[2x1 logical] 
10 

474625 
6.0849e-08 
[1x512 double] 


2.5 See Also 

interpl, griddedinterpolant, integraLg, funmin_g, meanMC_g, cubMC_g 


3 funmin_g 

l-D guaranteed global minimum value on [a,b] and the subset containing optimal solutions 

3.1 Syntax 

fmin = funmin_g(f) 

fmin = funmin_g(f,a,b,abstol,TolX) 

fmin = funmin_g(f,’a’,a,’b’,b,’abstor,abstol,’TolX’,TolX) 
fmin = funmin_g(f,in_param) 

[fmin, out_param] = funmin_g(f,...) 

3.2 Description 

fmin = funmin_g(f) finds minimum value of function f on the default interval [0,1] within the guaranteed 
absolute error tolerance of le-6 and the X tolerance of le-3. Default initial number of points is 100 and 
default cost budget is le7. Input f is a function handle. 

fmin = funmin_g(f,a,b,abstol,TolX) finds minimum value of function f with ordered input parameters that 
define the finite interval [a,b], a guaranteed absolute error tolerance abstol and a guaranteed X tolerance 
TolX. 

fmin = funmin_g(f,’a’, a, ’b’,b, ’abstol’,abstol,’TolX’,TolX) finds minimum value of function f on the interval 
[a,b] with a guaranteed absolute error tolerance abstol and a guaranteed X tolerance TolX. All five field-value 
pairs are optional and can be supplied in different order. 

fmin = funmin_g(f,in_param) finds minimum value of function f on the interval [in_param.a,in_param.b] 
with a guaranteed absolute error tolerance in_param.abstol and a guaranteed X tolerance in_param.TolX. If 
a field is not specified, the default value is used. 

[fmin, out_param] = funmin_g(f,...) returns minimum value fmin of function f and an output structure 
out_param. 

Input Arguments 

• f — input function 

• in_param.a — left end point of interval, default value is 0 

• in_param.b — right end point of interval, default value is 1 

• in_param.abstol — guaranteed absolute error tolerance, default value is le-6. 

• in_param.TolX — guaranteed X tolerance, default value is le-3. 

Optional Input Arguments 

• in_param.nlo — lower bound of initial number of points we used, default value is 10 

• in_param.nhi — upper bound of initial number of points we used, default value is 1000 
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• in_param.nmax — cost budget, default value is le7. 


Output Arguments 


• out_param.f — input function 

• out.param.a — left end point of interval 

• out.param.b — right end point of interval 

• out.param.abstol — guaranteed absolute error tolerance 

• out.param.TolX — guaranteed X tolerance 

• out.param.nlo — a lower bound of initial number of points we use 

• out_param.nhi — an upper bound of initial number of points we use 

• out_param.nmax — cost budget 

• out_param.ninit — initial number of points we use 

• out_param.tau — latest value of tau 

• out_param.npoints — number of points needed to reach the guaranteed absolute error tolerance or the 
guaranteed X tolerance 

• out_param.exitflag — the state of program when exiting 
0 Success 

1 Number of points used is greater than out_param.nmax 

• out_param.errest — estimation of the absolute error bound 

• out_param.volumeX — the volume of intervals containing the point (s) where the minimum occurs 

• out_param.tauchange — it is 1 if out_param.tau changes, otherwise it is 0 

• out_param.intervals — the intervals containing point(s) where the minimum occurs. Each column 
indicates one interval where the first row is the left point and the second row is the right point. 


3.3 Guarantee 

If the function to be minimized, / satisfies the cone condition 


< 


b-i 


f- 


fib) - fia) 


6-1 


then the fmin output by this algorithm is guaranteed to satisfy 

I min / — fmin I < abstol. 


or 

provided the flag exitflag = 0. 


volumeX < TolX, 
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3.4 Examples 

Example 1 

f=@(x) (x-0.3).~2+l; [fmin,out_param] = funmin_g(f) 

7o Minimize function (x-0.3)~2+l with default input parameter. 


fmin = 


1.0000 


out_param = 


f: 
a: 
b: 

abstol: 
TolX: 
nlo: 
nhi: 
nmax: 
ninit: 
tau: 
exitflag: 
npoints: 

errest: 
volumeX: 
tauchange: 
intervals: 


@(x)(x-0.3).'2+1 
0 
1 

l.OOOOe-06 

l.OOOOe-03 

10 

1000 

10000000 

100 

197 

0 

6337 

6.1554e-07 

0.0015 

0 

[2x1 double] 


Example 2 


f=@(x) (x-0.3).'2+1; 

[fmin,out_parain] = funmin_g(f,-2,2,le-7,le-4,10,10,1000000) 

°/o Minimize function (x-0.3)'2+1 on [-2,2] with error tolerance le-4, X 
7 tolerance le-2, cost budget 1000000, lower bound of initial number of 
7 points 10 and upper bound of initial number of points 10 


fmin = 


1.0000 


out_param = 

a: -2 

abstol: l.OOOOe-07 
b: 2 

f: @(x)(x-0.3).'2+1 
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nhi; 

10 

nlo; 

10 

nmax: 

1000000 

TolX: 

l.OOOOe-04 

ninit: 

10 

tau: 

17 

exitflag: 

0 

npoints: 

18433 

errest: 

9.5464e-08 

volumeX: 

5.4175e-04 

tauchange: 

0 

intervals: 

[2x1 double] 


Example 3 


clear in_parain; in_param.a = -13; in_parain.b = 8; 
in_parEiin. abstol = le-7; in_parain.TolX = le-4; 
in_parEim.nlo = 10; in_parain.nhi = 100; 
in_parEim. nmax = 10~6; 

[fmin, out_parain] = funmin_g(f , in_param) 

°/o Minimize function (x-0.3)"2+l on [-13,8] with error tolerance le-7, X 
"/o tolerance le-4, cost budget 1000000, lower bound of initial number of 
°/o points 10 and upper bound of initial number of points 100 


fmin = 


1 


out_param = 


a: 

-13 

abstol: 

l.OOOOe-07 

b: 

8 

f: 

@(x)(x-0.3).'2+l 

nhi: 

100 

nlo: 

10 

nmax: 

1000000 

TolX: 

l.OOOOe-04 

ninit: 

91 

tau: 

179 

exitflag: 

0 

npoints: 

368641 

errest: 

7.1014e-08 

volumeX: 

5.2445e-04 

tauchange: 

0 

intervals: 

[2x1 double] 
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Example 4 


f=@(x) (x-0.3).~2+l; 

[fmin,out_parain] = funmin_g(f , ’ a’,-2,’b',2,’nhi ’ , 100,’nlo', 10, . . . 

’nmax’,le6,’abstol’,le-4,'TolX’,le-2) 

°/o Minimize function (x-0.3)"2+l on [-2,2] with error tolerance le-4, X 
y, tolerance le-2, cost budget 1000000, lower bound of initial number of 
7o points 10 and upper bound of initial number of points 100 


fmin = 


1.0000 


out_param = 


a: 

abstol: 
b: 
f: 
nhi: 
nlo: 
nmax: 
TolX: 
ninit: 
tau: 
exitflag: 
npoints: 

errest: 
volumeX: 
tauchange: 
intervals: 


-2 

l.OOOOe-04 

2 

@(x)(x-0.3).'2+1 
100 
10 

1000000 

0.0100 

64 

125 

0 

2017 

6.2273e-05 

0.0146 

0 

[2x1 double] 


3.5 See Also 

fminbnd, funappx_g, integraLg 
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4 integral _g 

l-D guaranteed function integration using trapezoidal rule 

4.1 Syntax 

q = integral_g(f) 
q = integral_g(f,a,b,abstol) 
q = integral_g(f,’a’,a,’b’,b,’abstor,abstol) 
q = integral_g(f,in_param) 

[q, out_param] = integraLg(f,...) 

4.2 Description 

q = integral_g(f) computes q, the definite integral of function f on the interval [a,b] by trapezoidal rule with 
in a guaranteed absolute error of le-6. Default starting number of sample points taken is 100 and default 
cost budget is le7. Input f is a function handle. The function y = f(x) should accept a vector argument x 
and return a vector result y, the integrand evaluated at each element of x. 

q = integral_g(f,a,b,abstol) computes q, the definite integral of function f on the finite interval [a,b] by 
trapezoidal rule with the ordered input parameters, and guaranteed absolute error tolerance abstol. 

q = integral_g(f,’a’,a,’b’,b,’abstol’,abstol) computes q, the definite integral of function f on the finite interval 
[a,b] by trapezoidal rule within a guaranteed absolute error tolerance abstol. All four field-value pairs are 
optional and can be supplied. 

q = integral_g(f,in_param) computes q, the definite integral of function f by trapezoidal rule within a 
guaranteed absolute error in_param.abstol. If a field is not specified, the default value is used. 

[q, out_param] = integral_g(f,...) returns the approximated integration q and output structure out_param. 


Input Arguments 

• f — input function 

• in_param.a — left end of the integral, default value is 0 

• in_param.b — right end of the integral, default value is 1 

• in_param.abstol — guaranteed absolute error tolerance, default value is le-6 

Optional Input Arguments 

• in_param.nlo — lowest initial number of function values used, default value is 10 

• in_param.nhi — highest initial number of function values used, default value is 1000 

• in_param.nmax — cost budget (maximum number of function values), default value is le7 

• in_param.maxiter — max number of iterations, default value is 1000 
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Output Arguments 


• q — approximated integral 

• out_param.f — input function 

• out_param.a — low end of the integral 

• out_param.b — high end of the integral 

• out_param.abstol — guaranteed absolute error tolerance 

• out_param.nlo — lowest initial number of function values 

• out_param.nhi — highest initial number of function values 

• out_param.nmax — cost budget (maximum number of function values) 

• out_param.maxiter — max number of iterations 

• out_param.ninit — initial number of points we use, computed by nlo and nhi 

• out_param.tauchange — it is true if the cone constant has been changed, false otherwise. See [1] for 
details. If true, you may wish to change the input in_param.ninit to a larger number. 

• out_param.tauchange — it is true if the cone constant has been changed, false otherwise. See [1] for 
details. If true, you may wish to change the input in_param.ninit to a larger number. 

• out_param.iter — number of iterations 

• out_param.npoints — number of points we need to reach the guaranteed absolute error tolerance abstol. 

• out_param.errest — approximation error dehned as the differences between the true value and the 
approximated value of the integral. 

• out_param.nstar — final value of the parameter defining the cone of functions for which this algorithm 
is guaranteed; nstar = ninit-2 initially and is increased as necessary 

• out_param.exit — the state of program when exiting 
0 Success 

1 Number of points used is greater than out_param.nmax 

2 Number of iterations is greater than out_param.maxiter 


4.3 Guarantee 

If the function to be integrated, / satisfies the cone condition 

nstar 


iiriii < 


2(6 — a) 

then the q output by this algorithm is guaranteed to satisfy 

rb 


f{b) - fia) 


b — a 


/ fix)dx - q 

J a 


< abstol. 


provided the flag exceedbudget = 0. And the upper bound of the cost is 


nstar =i= (6 — a)^Ya,i'{f') 
2 X abstol 


+ 2 X nstar + 4. 


15 














4.4 Examples 

Example 1 

f = @(x) x.~2; [q, out_parEim] = integral_g(f) 

7„ Integrate function x with default input parameter to make the error less 
7» than le-7. 

q = 

0.3333 


out_param = 


f; 
a; 
b; 

abstol: 
nlo: 
nhi: 
nmax: 
maxiter: 
ninit: 
tau: 

exceedbudget: 
tauchange: 
iter: 

q: 

npoints: 
errest: 


@(x)x.~2 

0 

1 

l.OOOOe-06 

10 

1000 

10000000 

1000 

100 

197 

0 

0 

2 

0.3333 

3565 

9.9688e-07 


Example 2 


[q, out_parcun] = integral_g(@(x) exp(-x.~2),'a',1,’b’,2,... 

’nlo’,100,’nhi’,10000,’abstol’,le-5,’nmax’,le7) 

7„ Integrate function x~2 with starting number of points 52, cost budget 
7o 10000000 and error tolerance le-8 


0.1353 


out_param = 


a; 

abstol: 
b; 
f: 

maxiter: 


1 

l.OOOOe-05 

2 

@(x)exp(-x.~2) 
1000 
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nhi: 

10000 

nlo: 

100 

nmax: 

10000000 

ninit: 

1000 

tau: 

1997 

exceedbudget; 

0 

tauchange: 

0 

iter: 

2 

q: 

0.1353 

npoints; 

2998 

errest; 

7.3718e-06 


4.5 See Also 

integral, quad, funappx_g, meanMC_g, cubMC_g, funmin_g 
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5 meanMC_g 

Monte Carlo method to estimate the mean of a random variable 

5.1 Syntax 

tmu = meanMC_g(Yrand) 

tmu = meanMC_g(Yrand,abstol,reltol,alpha) 

tmu = meanMC_g(Yrand,’abstor,abstol,’reltor,reltol,’alpha’,alpha) 
[tmu, out_param] = meanMC_g(Yrand,in_param) 


5.2 Description 

tmu = meanMC_g(Yrand) estimates the mean, mu, of a random variable Y to within a specified generalized 
error tolerance, tolfun:=max(abstol,reltol*| mu j), i.e., | mu - tmu | <= tolfun with probability at least 1- 
alpha, where abstol is the absolute error tolerance, and reltol is the relative error tolerance. Usually the reltol 
determines the accuracy of the estimation, however, if the j mu j is rather small, the abstol determines the 
accuracy of the estimation. The default values are abstol=le-2, reltol=le-l, and alpha=l%. Input Yrand is 
a function handle that accepts a positive integer input n and returns an n x 1 vector of IID instances of the 
random variable Y. 

tmu = meanMC_g(Yrand,abstol,reltol,alpha) estimates the mean of a random variable Y to within a 
specified generalized error tolerance tolfun with guaranteed confidence level 1-alpha using all ordered parsing 
inputs abstol, reltol, alpha. 

tmu = meanMC_g(Yrand,’abstol’,abstol,’reltol’,reltol,’alpha’,alpha) estimates the mean of a random vari¬ 
able Y to within a specified generalized error tolerance tolfun with guaranteed confidence level 1-alpha. All 
the field-value pairs are optional and can be supplied in different order, if a field is not supplied, the default 
value is used. 

[tmu, out_param] = meanMC_g(Yrand,in_param) estimates the mean of a random variable Y to within a 
specified generalized error tolerance tolfun with the given parameters in_param and produce the estimated 
mean tmu and output parameters out_param. If a field is not specified, the default value is used. 

Input Arguments 

• Yrand — the function for generating n IID instances of a random variable Y whose mean we want to 
estimate. Y is often defined as a function of some random variable X with a simple distribution. The 
input of Yrand should be the number of random variables n, the output of Yrand should be n function 
values. For example, if Y = X."2 where X is a standard uniform random variable, then one may define 
Yrand = @(n) rand(n,l)."2. 

• in_param.abstol — the absolute error tolerance, which should be positive, default value is le-2. 

• in_param.reltol — the relative error tolerance, which should be between 0 and 1, default value is le-1. 

• in_param.alpha — the uncertainty, which should be a small positive percentage, default value is 1%. 

Optional Input Arguments 

• in_param.fudge — standard deviation inflation factor, which should be larger than 1, default value is 

1 . 2 . 
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• in_param.nSig — initial sample size for estimating the sample variance, which should be a moderate 
large integer at least 30, the default value is led. 

• in_param.nl — initial sample size for estimating the sample mean, which should be a moderate large 
positive integer at least 30, the default value is led. 

• in_param.tbudget — the time budget in seconds to do the two-stage estimation, which should be 
positive, the default value is 100 seconds. 

• in_param.nbudget — the sample budget to do the two-stage estimation, which should be a large positive 
integer, the default value is le9. 

Output Arguments 


• tmu — the estimated mean of Y. 

• out_param.tau — the iteration step. 

• out_param.n — the sample size used in each iteration. 

• out_param.nremain — the remaining sample budget to estimate mu. It was calculated by the sample 
left and time left. 

• out_param.ntot — total sample used. 

• out_param.hmu — estimated mean in each iteration. 

• out_param.tol — the reliable upper bound on error for each iteration. 

• out_param.var — the sample variance. 

• out_param.exit — the state of program when exiting. 

0 Success 

1 Not enough samples to estimate the mean 

• out_param.kurtmax — the upper bound on modified kurtosis. 

• out_param.time — the time elapsed in seconds. 

• out_param.flag — parameter checking status 
1 checked by meanMC_g 

5.3 Guarantee 

This algorithm attempts to calculate the mean, mu, of a random variable to a prescribed error tolerance, 
tolfun:= max(abstol,reltol*| mu |), with guaranteed confidence level 1-alpha. If the algorithm terminated 
without showing any warning messages and provide an answer tmu, then the follow inequality would be 
satisfied: Pr(| mu - tmu | <= tolfun) >= 1-alpha The cost of the algorithm, N_tot, is also bounded above 

by N_up, which is defined in terms of abstol, reltol, nSig, nl, fudge, kurtmax, beta. And the following 
inequality holds: Pr (N_tot <= N_up) >= 1-beta Please refer to our paper for detailed arguments and 
proofs. 
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5.4 Examples 

Example 1 


7o Calculate the mean of x~2 when x is uniformly distributed in 

7o [0 1], with the absolute error tolerance = le-3 and uncertainty 57,. 

in_param.reltol=0; in_param.abstol = le-3; in_param.reltol = 0; 
in_param.alpha = 0.05; Yrand=@(n) rand(n,1).~2; 
tmu = meanMC_g(Yrand,in_parain) 


tmu = 


0.3331 


Example 2 

7o Calculate the mean of exp(x) when x is uniformly distributed in 
7« [0 1], with the absolute error tolerance le-3. 

tmu = meanMC_g(@(n)exp(rand(n,1)),le-3,0) 

tmu = 

1.7185 


Example 3 


7o Calculate the mean of cos(x) when x is uniformly distributed in 
7o [0 1], with the relative error tolerance le-2 and uncertainty 0.05. 

tmu = meanMC_g(@(n)cos(rand(n,1)),’reltol’,le-2abstol’,0,... 
’alpha’,0.05) 


tmu = 


0.8415 


5.5 See Also 

funappx_g, integraLg, cubMC_g, meanMCBer_g, cubSoboLg, cubLattice_g 
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6 meanMCBer_g 

Monte Carlo method to estimate the mean of a Bernoulli random variable to within a specified absolute 
error tolerance with guaranteed confidence level 1-alpha. 


6.1 Syntax 

pHat = meanMCBer_g(Yrand) 

pHat = meanMCBer_g(Yrand,abstol,alpha,nmax) 

pHat = meanMCBer_g(Yrand,’abstor,abstol,’alpha’,alpha,’nmax’,nmax) 
[pHat, out_param] = meanMCBer_g(Yrand,in_param) 


6.2 Description 

pHat = meanMCBer_g(Yrand) estimates the mean of a Bernoulli random variable Y to within a specified 
absolute error tolerance with guaranteed confidence level 99%. Input Yrand is a function handle that accepts 
a positive integer input n and returns a n x 1 vector of HD instances of the Bernoulli random variable Y. 

pHat = meanMCBer_g(Yrand,abstol,alpha,nmax) estimates the mean of a Bernoulli random variable Y to 
within a specified absolute error tolerance with guaranteed confidence level 1-alpha using all ordered parsing 
inputs abstol, alpha and nmax. 

pHat = meanMCBer_g(Yrand,’abstol’,abstol,’alpha’,alpha,’nmax’,nmax) estimates the mean of a Bernoulli 
random variable Y to within a specified absolute error tolerance with guaranteed conhdence level 1-alpha. 
All the field-value pairs are optional and can be supplied in different order. 

[pHat, out_param] = meanMCBer_g(Yrand,in_param) estimates the mean of a Bernoulli random variable Y 
to within a specified absolute error tolerance with the given parameters in_param and produce the estimated 
mean pHat and output parameters out_param. 

Input Arguments 


• Yrand — the function for generating HD instances of a Bernoulli random variable Y whose mean we 
want to estimate. 

• pHat — the estimated mean of Y. 

• in_param.abstol — the absolute error tolerance, the default value is le-2. 

• in_param.alpha — the uncertainty, the default value is 1%. 

• in_param.nmax — the sample budget, the default value is le9. 

Output Arguments 

• out_param.n — the total sample used. 

• out_param.time — the time elapsed in seconds. 

• out_param.exit — the state of program when exiting. 

0 Success 

1 Not enough samples to estimate p with guarantee 
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6.3 Guarantee 


If the sample size is calculated according Hoeffding’s inequality, which equals to 
ceil(log(2/out_param.alpha)/(2*out_param.abstor2)), then the following inequality must be satisfied: 

Pr(| p - pHat I <= abstol) >= 1-alpha. Here p is the true mean of Yrand, and pHat is the output of 

MEANMCBER_G. Also, the cost is deterministic. 

6.4 Examples 

Example 1 


7o Calculate the mean of a Bernoulli random variable with true p=l/90, 

7„ absolute error tolerance le-3 and uncertainty 0.01. 

in_parcun.abstol = le-3; in_param.alpha = 0.01; in_param.nmax = le9; 
p=l/9; Yrand=@(n) rand(n,l)<p; 
pHat = meanMCBer_g(Yrand,in_param) 


pHat = 


0.1113 


Example 2 


7o Using the same function as exaimple 1, with the absolute error tolerance 
1 le-4. 


pHat = meanMCBer_g(Yrand,le-4) 


pHat = 


0.1111 


Example 3 


7o Using the same function as exaimple 1, with the absolute error tolerance 
7o le-2 and uncertainty 0.05. 

pHat = meanMCBer_g(Yrand,’abstol’,le-2,’alpha’,0.05) 


pHat = 


0.1118 


6.5 See Also 

funappx_g, integraLg, cubMC_g, meanMC_g, cubLattice_g, cubSoboLg 
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7 cubMC_g 

Monte Carlo method to evaluate a multidimensional integral 


7.1 Syntax 

[Q,out_param] = cubMC_g(f,hyperbox) 

Q = cubMC_g(f,hyperbox,measure,abstol,reltol,alpha) 

Q = cubMC_g(f,hyperbox, ’measure’,measure, ’abstor,abstol, ’reltol’ ,reltol, ’alpha’,alpha) 

[Q out_param] = cubMC_g(f,hyperbox,in_param) 

7.2 Description 

[Q,out_param] = cubMC_g(f,hyperbox) estimates the integral of f over hyperbox to within a specified 
generalized error tolerance, tolfun = max(abstol, reltol*| I |), i.e., | I - Q | <= tolfun with probability at 
least 1-alpha, where abstol is the absolute error tolerance, and reltol is the relative error tolerance. Usually 
the reltol determines the accuracy of the estimation, however, if the | I | is rather small, the abstol determines 
the accuracy of the estimation. The default values are abstol=le-2, reltol=le-l, and alpha=l%. Input f is 
a function handle that accepts an n x d matrix input, where d is the dimension of the hyperbox, and n is 
the number of points being evaluated simultaneously. The input hyperbox is a 2 x d matrix, where the first 
row corresponds to the lower limits and the second row corresponds to the upper limits. 

Q = cubMC_g(f,hyperbox,measure,abstol,reltol,alpha) estimates the integral of function f over hyperbox 
to within a specified generalized error tolerance tolfun with guaranteed confidence level 1-alpha using all 
ordered parsing inputs f, hyperbox, measure, abstol, reltol, alpha, fudge, nSig, nl, tbudget, nbudget, flag. 
The input f and hyperbox are required and others are optional. 

Q = cubMC_g(f,hyperbox,’measure’,measure,’abstol’,abstol,’reltol’,reltol,’alpha’,alpha) estimates the inte¬ 
gral of f over hyperbox to within a specified generalized error tolerance tolfun with guaranteed confidence 
level 1-alpha. All the field-value pairs are optional and can be supplied in different order. If an input is not 
specified, the default value is used. 

[Q out_param] = cubMC_g(f,hyperbox,in_param) estimates the integral of f over hyperbox to within a spec¬ 
ified generalized error tolerance tolfun with the given parameters in_param and produce output parameters 
out_param and the integral Q. 

Input Arguments 


• f — the integrand. 

• hyperbox — the integration hyperbox. The default value is [zeros(l,d); ones(l,d)], the default d is 1. 

• in_param.measure — the measure for generating the random variable, the default is ’uniform’. The 
other measure could be handled is ’normal’/’Gaussian’. The input should be a string type, hence with 
quotes. 

• in_param.abstol — the absolute error tolerance, the default value is le-2. 

• in_param.reltol — the relative error tolerance, the default value is le-1. 

• in_param.alpha — the uncertainty, the default value is 1%. 
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Optional Input Arguments 


• in_param.fudge — the standard deviation inflation factor, the default value is 1.2. 

• in_param.nSig — initial sample size for estimating the sample variance, which should be a moderate 
large integer at least 30, the default value is led. 

• in_param.nl — initial sample size for estimating the sample mean, which should be a moderate large 
positive integer at least 30, the default value is led. 

• in_param.tbudget — the time budget to do the estimation, the default value is 100 seconds. 

• in_param.nbudget — the sample budget to do the estimation, the default value is le9. 

• in_param.flag — the value corresponds to parameter checking status. 

0 not checked 

1 checked by meanMC_g 

2 checked by cubMC.g 

Output Arguments 


• Q — the estimated value of the integral. 

• out_param.n — the sample size used in each iteration. 

• out_param.ntot — total sample used. 

• out_param.nremain — the remaining sample budget to estimate 1. It was calculated by the sample left 
and time left. 

• out_param.tau — the iteration step. 

• out_param.hmu — estimated integral in each iteration. 

• out_param.tol — the reliable upper bound on error for each iteration. 

• out_param.kurtmax — the upper bound on modified kurtosis. 

• out_param.time — the time elapsed in seconds. 

• out_param.var — the sample variance. 

• out_param.exit — the state of program when exiting. 

0 success 

I Not enough samples to estimate the mean 
10 hyperbox does not contain numbers 

II hyperbox is not 2 x d 

12 hyperbox is only a point in one direction 

13 hyperbox is infinite when measure is ’uniform’ 

14 hyperbox is not doubly infinite when measure is ’normal’ 
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7.3 Guarantee 


This algorithm attempts to calculate the integral of function f over a hyperbox to a prescribed error tolerance 
tolfun:= max(abstol,reltol*| I |) with guaranteed confidence level 1-alpha. If the algorithm terminated 
without showing any warning messages and provide an answer Q, then the follow inequality would be 
satisfied: 

Pr(| Q - I I <= tolfun) >= 1-alpha 

The cost of the algorithm, N_tot, is also bounded above by N_up, which is a function in terms of abstol, 
reltol, nSig, nl, fudge, kurtmax, beta. And the following inequality holds: 

Pr (N_tot <= N_up) >= 1-beta 

Please refer to our paper for detailed arguments and proofs. 


7.4 Examples 

Example 1 

"/o Estimate the integral with integrand f(x) = sin(x) over the interval 

y. [1;2] 

f = @(x) sin(x); interval = [1;2]; 

Q = cubMC_g(f,interval,’uniform’,le-3,le-2) 

q = 


0.9564 


Example 2 

7o Estimate the integral with integrand f(x) = exp(-xl~2-x2~2) over the 
y hyperbox [0 0;1 1], where x is a vector x = [xl x2] . 

f = @(x) exp(-x(:,1)."2-x(:,2).~2); hyperbox = [0 0;1 1]; 

Q = cubMC_g(f,hyperbox,’measure’,’uniform’,’abstol’,le-3,... 

’reltol’,le-13) 


Q = 


0.5574 


Example 3 

y Estimate the integral with integrand f(x) = 2~d*prod(xl*x2*...*xd) + 
y 0.555 over the hyperbox [zeros(l,d);ones(l,d)], where x is a vector x = 
y [xl x2... xd]. 


d = 3;f = @(x) 2~d*prod(x,2)+0.555; hyperbox = [zeros(l,d);ones(l,d)]; 
in_param.abstol = le-3; in_param.reltol=le-3; 

Q = cubMC_g(f,hyperbox,in_param) 
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1.5549 


Q = 


Example 4 


7, Estimate the integral with integrand f(x) = exp(-xl~2-x2~2) in the 
7o hyperbox [-inf -inf; inf inf] , where x is a vector x = [xl x2] . 

f = @(x) exp(-x(:,1)."2-x(:,2).~2); hyperbox = [-inf -inf;inf inf]; 
Q = cubMC_g(f,hyperbox,’normal’,0,le-2) 

Q = 


0.3328 


7.5 See Also 

funappx_g, integraLg, meanMC_g, meanMCBer_g, cubLattice_g, cubSoboLg 
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8 cubLattice_g 

Quasi-Monte Carlo method using rank-1 Lattices cubature over a d-dimensional region to integrate within 
a specified generalized error tolerance with guarantees under Fourier coefficients cone decay assumptions. 


8.1 Syntax 

[q,out_param] = cubLattice_g(f,hyperbox) 
q = cubLattice_g(f,hyperbox,measure,abstol,reltol) 

q = cubLattice_g(f,hyperbox, ’measure’,measure, ’abstol’ ,abstol, ’reltol’ ,reltol) 
q = cubLattice_g(f,hyperbox,in_param) 


8.2 Description 

[q,out_param] = cubLattice_g(f,hyperbox) estimates the integral of f over the d-dimensional region de¬ 
scribed by hyperbox, and with an error guaranteed not to be greater than a specific generalized error 
tolerance, tolfun:=max(abstol,reltol*| integral(f) |). Input f is a function handle, f should accept an n x d 
matrix input, where d is the dimension and n is the number of points being evaluated simultaneously. The 
input hyperbox is a 2 x d matrix, where the first row corresponds to the lower limits and the second row 
corresponds to the upper limits of the integral. Given the construction of our Lattices, d must be a positive 
integer with l<=d<=250. 

q = cubLattice_g(f, hyperbox,measure,abstol,reltol) estimates the integral of f over the hyperbox. The 
answer is given within the generalized error tolerance tolfun. All parameters should be input in the order 
specified above. If an input is not specified, the default value is used. Note that if an input is not specified, 
the remaining tail cannot be specified either. Inputs f and hyperbox are required. The other optional inputs 
are in the correct order: measure,abstol,reltol,shift,mmin,mmax,fudge,transform,toltype and theta. 

q = cubLattice_g(f, hyperbox,’measure’,measure,’abstol’,abstol,’reltol’,reltol) estimates the integral of f over 
the hyperbox. The answer is given within the generalized error tolerance tolfun. All the field-value pairs are 
optional and can be supplied in any order. If an input is not specified, the default value is used. 

q = cubLattice_g (f,hyperbox,in_param) estimates the integral of f over the hyperbox. The answer is given 
within the generalized error tolerance tolfun. 

Input Arguments 

• f — the integrand whose input should be a matrix n x d where n is the number of data points and d 
the dimension, which cannot be greater than 250. By default f is f=@ x."2. 

• hyperbox — the integration region defined by its bounds. It must be a 2 x d matrix, where the first 
row corresponds to the lower limits and the second row corresponds to the upper limits of the integral. 
The default value is [0;I]. 

• in_param.measure — for f(x)*mu(dx), we can define mu(dx) to be the measure of a uniformly dis¬ 
tributed random variable in they hyperbox or normally distributed with covariance matrix I_d. The 
only possible values are ’uniform’ or ’normal’. For ’uniform’, the hyperbox must be a finite volume 
while for ’normal’, the hyperbox can only be defined as (-Inf,Inf)"d. By default it is ’uniform’. 

• in_param.abstol — the absolute error tolerance, abstol>=0. By default it is Ie-4. 

• in_param.reltol — the relative error tolerance, which should be in [0,1]. Default value is le-2. 
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Optional Input Arguments 


• in_param.shift — the Rank-1 lattices can be shifted to avoid the origin or other particular points. By 
default we consider a uniformly [0,1) random shift. 

• in_param.mmin — the minimum number of points to start is 2'mmin. The cone condition on the 
Fourier coefficients decay requires a minimum number of points to start. The advice is to consider at 
least mmin=10. mmin needs to be a positive integer with mmin<=mmax. By default it is 10. 

• in_param.mmax — the maximum budget is 2''mmax. By construction of our Lattices generator, mmax 
is a positive integer such that mmin<=mmax<=26. The default value is 24. 

• in_param.fudge — the positive function multiplying the finite sum of Fast Fourier coefficients specified 
in the cone of functions. This input is a function handle. The fudge should accept an array of nonneg¬ 
ative integers being evaluated simultaneously. For more technical information about this parameter, 
refer to the references. By default it is @(m) 5*2."-m. 

• in_param.transform — the algorithm is defined for continuous periodic functions. If the input function 
f is not, there are 5 types of transform to periodize it without modifying the result. By default it is 
the Baker’s transform. The options are: 

’id’ : no transformation. 

’Baker’ : Baker’s transform or tent map in each coordinate. Preserving only continuity but simple to 
compute. Chosen by default. 

’CO’ : polynomial transformation only preserving continuity. 

’Cl’ : polynomial transformation preserving the first derivative. 

’Clsin’ : Sidi’s transform with sine, preserving the first derivative. This is in general a better option 
than ’Cl’. 

• in_param.toltype — this is the generalized tolerance function. There are two choices, ’max’ which 
takes max(abstol,reltol*| integral(f) | ) and ’comb’ which is the linear combination theta*abstol-|-(l- 
theta)*reltol*| integral(f) | . Theta is another parameter to be specified with ’comb’(see below). For 
pure absolute error, either choose ’max’ and set reltol = 0 or choose ’comb’ and set theta = 1. For 
pure relative error, either choose ’max’ and set abstol = 0 or choose ’comb’ and set theta = 0. Note 
that with ’max’, the user can not input abstol = reltol = 0 and with ’comb’, if theta = 1 abstol con 
not be 0 while if theta = 0, reltol can not be 0. By default toltype is ’max’. 

• in_param.theta — this input is parametrizing the toltype ’comb’. Thus, it is only active when the 
toltype chosen is ’comb’. It establishes the linear combination weight between the absolute and relative 
tolerances theta*abstol-|-(1-theta)*reltol*| integral(f) |. Note that for theta = 1, we have pure absolute 
tolerance while for theta = 0, we have pure relative tolerance. By default, theta=l. 

Output Arguments 


• q — the estimated value of the integral. 

• out_param.d — dimension over which the algorithm integrated. 

• out_param.n — number of Rank-1 lattice points used for computing the integral of f. 

• out_param.bound_err — predicted bound on the error based on the cone condition. If the function lies 
in the cone, the real error will be smaller than generalized tolerance. 

• out_param.time — time elapsed in seconds when calling cubLattice_g. 
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• out_param.exitfiag — this is a binary vector stating whether warning flags arise. These flags teii about 
which conditions make the flnai resnit certainiy not guaranteed. One flag is considered arisen when its 
vaiue is 1. The foiiowing iist expiains the flags in the respective vector order: 

1 If reaching overbudget. It states whether the max budget is attained without reaching the guar¬ 
anteed error tolerance. 

2 If the function lies outside the cone. In this case, results are not guaranteed. Note that this 
parameter is computed on the transformed function, not the input function. For more information 
on the transforms, check the input parameter in_param.transform; for information about the cone 
definition, check the article mentioned below. 

8.3 Guarantee 

This algorithm computes the integral of real valued functions in dimension d with a prescribed generalized 
error tolerance. The Fourier coefficients of the integrand are assumed to be absolutely convergent. If the 
algorithm terminates without warning messages, the output is given with guarantees under the assumption 
that the integrand lies inside a cone of functions. The guarantee is based on the decay rate of the Fourier 
coefficients. For more details on how the cone is defined, please refer to the references below. 


8.4 Examples 

Example 1 

"/o Estimate the integral with integrand f (x) = xl.*x2 in the interval 

7 . [ 0 , 1 )- 2 : 

f = @(x) prod(x,2); hyperbox = [zeros(l,2);ones(l,2)]; 
q = cubLattice_g(f,hyperbox,'uniform',le-5,0,'transform','Clsin') 

q = 


0.2500 


Example 2 

7, Estimate the integral with integrand f(x) = xl. ~2 . *x2 . ~2 . *x3. "2 
7o in the interval R~3 where xl, x2 and x3 are normally distributed: 

f = @(x) x(:,1).~2.*x(:,2).~2.*x(:,3).~2; hyperbox = [-inf(1,3);inf(1,3)]; 
q = cubLattice_g(f,hyperbox,'normal',le-3,le-3,'transform','Clsin') 

q = 


1.0000 


Example 3 

7o Estimate the integral with integrand f(x) = exp(-xl~2-x2~2) in the 
7o interval [-1,2)~2: 

f = @(x) exp(-x(:,1).~2-x(:,2)."2); hyperbox = [-ones(l,2);2*ones(l,2)]; 
q = cubLattice_g(f,hyperbox,'uniform',le-3,le-2,'transform','Cl') 
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2.6532 


q = 


Example 4 


7, Estimate the price of an European call with 30=100, K=100, r=sigma~2/2, 

7o sigma=0.05 and T=l. 

f = @(x) exp(-0.05~2/2)*max(100*exp(0.05*x)-100,0); hyperbox = [-inf(1,1);inf(1,1)] 
q = cubLattice_g(f,hyperbox,'normal’,le-4,le-2,'transform','Clsin') 

q = 


2.0563 


Example 5 


7, Estimate the integral with integrand f(x) = 8*xl.*x2.*x3.*x4.*x5 in the 
7o interval [0,1) ~5 with pure absolute error le-5. 

f = @(x) 8*prod(x,2); hyperbox = [zeros(l,5);ones(l,5)]; 
q = cubLattice_g(f,hyperbox,'uniform',le-5,0) 

q = 


0.2500 
Example 6 


7o Estimate the integral with integrand f(x) = 3./(5-4*(cos(2*pi*x))) in the interval 
7o [0,1) with pure absolute error le-5. 

f = @(x) 3./(5-4*(cos(2*pi*x))); hyperbox = [0;1]; 
q = cubLattice_g(f,hyperbox,'uniform',le-5,0,'transform','id') 

q = 


1.0000 

8.5 See Also 

cubSoboLg, cubMC_g, meanMC_g, meanMCBer_g, integraLg 
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9 cubSoboLg 

Quasi-Monte Carlo method using Sobol’ cubature over the d-dimensional region to integrate within a specified 
generalized error tolerance with guarantees under Walsh-Fourier coefficients cone decay assumptions 


9.1 Syntax 

[q,out_param] = cubSobol_g(f,hyperbox) 
q = cubSobol_g(f,hyperbox,measure,abstol,reltol) 

q = cubSobol_g(f,hyperbox,’measure’,measure,’abstor,abstol,’reltor,reltol) 
q = cubSobol_g(f,hyperbox,in_param) 


9.2 Description 

[q,out_param] = cubSoboLg(f, hyperbox) estimates the integral of f over the d-dimensional region described 
by hyperbox, and with an error guaranteed not to be greater than a specific generalized error tolerance, 
tolfun:=max(abstol,reltol*| integral(f) |). Input f is a function handle, f should accept an n x d matrix 
input, where d is the dimension and n is the number of points being evaluated simultaneously. The input 
hyperbox is a 2 x d matrix, where the first row corresponds to the lower limits and the second row corresponds 
to the upper limits of the integral. Given the construction of Sobol’ sequences, d must be a positive integer 
with l<=d<=llll. 

q = cubSoboLg(f,hyperbox,measure,abstol,reltol) estimates the integral of f over the hyperbox. The answer 
is given within the generalized error tolerance tolfun. All parameters should be input in the order specified 
above. If an input is not specified, the default value is used. Note that if an input is not specified, the 
remaining tail cannot be specified either. Inputs f and hyperbox are required. The other optional inputs are 
in the correct order: measure,abstol,reltol,mmin,mmax,fudge,toltype and theta. 

q = cubSoboLg(f,hyperbox,’measure’,measure,’abstol’,abstol,’reltol’,reltol) estimates the integral of f over 
the hyperbox. The answer is given within the generalized error tolerance tolfun. All the field-value pairs are 
optional and can be supplied in any order. If an input is not specified, the default value is used. 

q = cubSobol_g(f,hyperbox,in_param) estimates the integral of f over the hyperbox. The answer is given 
within the generalized error tolerance tolfun. 

Input Arguments 

• f — the integrand whose input should be a matrix n x d where n is the number of data points and d 
the dimension, which cannot be greater than 1111. By default f is f=@ x."2. 

• hyperbox — the integration region defined by its bounds. It must be a 2 x d matrix, where the first 
row corresponds to the lower limits and the second row corresponds to the upper limits of the integral. 
The default value is [0;1]. 

• in_param.measure — for f(x)*mu(dx), we can define mu(dx) to be the measure of a uniformly dis¬ 
tributed random variable in the hyperbox or normally distributed with covariance matrix I_d. The 
only possible values are ’uniform’ or ’normal’. For ’uniform’, the hyperbox must be a finite volume 
while for ’normal’, the hyperbox can only be defined as (-Inf,Inf)"d. By default it is ’uniform’. 

• in_param.abstol — the absolute error tolerance, abstol>=0. By default it is le-4. 

• in_param.reltol — the relative error tolerance, which should be in [0,1]. Default value is le-2. 
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Optional Input Arguments 


• iii_param.mmiii — the minimum number of points to start is 2'mmin. The cone condition on the 
Fourier coefficients decay requires a minimum number of points to start. The advice is to consider at 
least mmin=10. mmin needs to be a positive integer with mmin<=mmax. By default it is 10. 

• in_param.mmax — the maximum budget is 2'mmax. By construction of the Sobol’ generator, mmax 
is a positive integer such that mmin<=mmax<=53. The default value is 24. 

• in_param.fudge — the positive function multiplying the finite sum of Fast Walsh Fourier coefficients 
specified in the cone of functions. This input is a function handle. The fudge should accept an array 
of nonnegative integers being evaluated simultaneously. For more technical information about this 
parameter, refer to the references. By default it is @(m) 5*2.*-m. 

• in_param.toltype — this is the generalized tolerance function. There are two choices, ’max’ which 
takes max(abstol,reltol*| integral(f) | ) and ’comb’ which is the linear combination theta*abstol+(l- 
theta)*reltol*| integral(f) | . Theta is another parameter to be specified with ’comb’(see below). For 
pure absolute error, either choose ’max’ and set reltol = 0 or choose ’comb’ and set theta = 1. For 
pure relative error, either choose ’max’ and set abstol = 0 or choose ’comb’ and set theta = 0. Note 
that with ’max’, the user can not input abstol = reltol = 0 and with ’comb’, if theta = 1 abstol con 
not be 0 while if theta = 0, reltol can not be 0. By default toltype is ’max’. 

• in_param.theta — this input is parametrizing the toltype ’comb’. Thus, it is only active when the 
toltype chosen is ’comb’. It establishes the linear combination weight between the absolute and relative 
tolerances theta*abstol+(l-theta)*reltol*| integral(f) |. Note that for theta = 1, we have pure absolute 
tolerance while for theta = 0, we have pure relative tolerance. By default, theta=l. 

Output Arguments 

• q — the estimated value of the integral. 

• out_param.d — dimension over which the algorithm integrated. 

• out_param.n — number of Sobol’ points used for computing the integral of f. 

• out_param.bound_err — predicted bound on the error based on the cone condition. If the function lies 
in the cone, the real error will be smaller than generalized tolerance. 

• out_param.time — time elapsed in seconds when calling cubSoboLg. 

• out_param.exitflag — this is a binary vector stating whether warning flags arise. These flags tell about 
which conditions make the final result certainly not guaranteed. One flag is considered arisen when its 
value is 1. The following list explains the flags in the respective vector order: 

1 If reaching overbudget. It states whether the max budget is attained without reaching the guar¬ 
anteed error tolerance. 

2 If the function lies outside the cone. In this case, results are not guaranteed. For more information 
about the cone definition, check the article mentioned below. 

9.3 Guarantee 

This algorithm computes the integral of real valued functions in dimension d with a prescribed generalized 
error tolerance. The Walsh-Fourier coefficients of the integrand are assumed to be absolutely convergent. 
If the algorithm terminates without warning messages, the output is given with guarantees under the as¬ 
sumption that the integrand lies inside a cone of functions. The guarantee is based on the decay rate of the 
Walsh-Fourier coefficients. For more details on how the cone is defined, please refer to the references below. 
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9.4 Examples 

Example 1 

7„ Estimate the integral with integrand f(x) = xl.*x2 in the interval 

1 [ 0 , 1 )‘ 2 : 

f = @(x) prod(x,2); hyperbox = [zeros(l,2);ones(l,2)]; 
q = cubSobol_g(f,hyperbox,’uniform’,le-5,0) 

q = 


0.2500 


Example 2 

7„ Estimate the integral with integrand f (x) = xl. ~2 . *x2 . ~2 . *x3. ”2 
7„ in the interval R~3 where xl, x2 and x3 are normally distributed: 

f = @(x) x(:,1).~2.*x(:,2).~2.*x(:,3).~2; hyperbox = [-inf(1,3);inf(1,3)]; 
q = cubSobol_g(f,hyperbox,’normal’,le-3,le-3) 

q = 


1.0004 


Example 3 

7, Estimate the integral with integrand f(x) = exp(-xl~2-x2~2) in the 
7o interval [-1,2)~2: 

f = @(x) exp(-x(:,1).~2-x(:,2)."2); hyperbox = [-ones(l,2);2*ones(l,2)]; 
q = cubSobol_g(f,hyperbox,’uniform’,le-3,le-2) 

q = 


2.6532 


Example 4 

7o Estimate the price of an European call with 30=100, K=100, r=sigma~2/2, 

7o sigma=0.05 and T=l. 

f = @(x) exp(-0.05~2/2)*max(100*exp(0.05*x)-100,0); hyperbox = [-inf(1,1);inf(1,1)] 
q = cubSobol_g(f,hyperbox,’normal’,le-4,le-2) 

q = 


2.0552 
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Example 5 


7„ Estimate the integral with integrand f(x) = 8*xl.*x2.*x3.*x4.*x5 in the 
7„ interval [0,1) ~5 with pure absolute error le-5. 

f = @(x) 8*prod(x,2); hyperbox = [zeros(l,5);ones(l,5)]; 
q = cubSobol_g(f,hyperbox,'uniform’,le-5,0) 

q = 


0.2500 


9.5 See Also 

cubLattice_g, cubMC_g, meanMC_g, meanMCBer_g, integraLg 
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